

#include "ReconstructionSystem.h"
#include "FEProblem.h"
#include "VFVProblem.h"
#include "Factory.h"
#include "AuxKernel.h"
#include "AuxScalarKernel.h"
#include "MaterialData.h"
#include "Assembly.h"
#include "GeometricSearchData.h"
#include "ComputeNodalAuxVarsThread.h"
#include "ComputeNodalAuxBcsThread.h"
#include "ComputeElemAuxVarsThread.h"
#include "ComputeElemAuxBcsThread.h"
#include "Parser.h"
#include "TimeIntegrator.h"
#include "Conversion.h"

#include "NonlinearSystem.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/node_range.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/default_coupling.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/fe_interface.h"
#include "libmesh/enum_solver_type.h"

#include <chrono>
using namespace std::chrono;

void assemble_reconstrunction(EquationSystems & es, const std::string & system_name)
{
    // 获取EquationSystem的引用

    auto start = steady_clock::now();

    LinearImplicitSystem & _sys = es.get_system<LinearImplicitSystem>(system_name);

    // 从EquationSystem中取出MeshBase
    const auto & mesh = es.get_mesh();

    // Mesh的维数，用来创建FEBase
    const unsigned int dim = mesh.mesh_dimension();

    // 从LinearImplictSystem中取出DofMap
    const auto& dof_map = _sys.get_dof_map();

    // 取出FEType，用来创建FEBase
    auto fe_type = _sys.variable_type(0);

    // 单元的FEBase，存储和更新积分的基函数
    std::unique_ptr<FEBase> fe(FEBase::build(dim, fe_type));
    // 面上的FEBase，存储和更新左单元在面上积分点的基函数
    std::unique_ptr<FEBase> fe_elem_face(FEBase::build(dim, fe_type));
    // 面上的FEBase，存储和更新右单元在面上积分点的基函数
    std::unique_ptr<FEBase> fe_neighbor_face(FEBase::build(dim, fe_type));

    // 创建Cell的Gauss积分，默认精度2p+1
    // TODO: 积分精度只需要p+1
    QGauss qrule (dim, fe_type.default_quadrature_order());
    fe->attach_quadrature_rule (&qrule);

    // 创建Face的Gauss积分，维度=dim-1，默认精度2p+1
    QGauss qface (dim-1, fe_type.default_quadrature_order());
    fe_elem_face->attach_quadrature_rule(&qface);
    fe_neighbor_face->attach_quadrature_rule(&qface);

    // 从Cell的fe中取出Gauss积分点上的值，如权重，基函数，基函数的梯度
    // 注意：这里要用引用！
    const auto& JxW = fe->get_JxW();
    const auto& phi = fe->get_phi();
    const auto& dphi = fe->get_dphi();

    // Face的左值
    const auto&  phi_face = fe_elem_face->get_phi();
    const auto& dphi_face = fe_elem_face->get_dphi();

    // Face的右值
    const auto&  phi_neighbor_face = fe_neighbor_face->get_phi();
    const auto& dphi_neighbor_face = fe_neighbor_face->get_dphi();

    // Face上积分权重    
    const auto& JxW_face = fe_elem_face->get_JxW();

    // 法向量
    const auto& qface_normals = fe_elem_face->get_normals();

    // Face上的Gauss积分点坐标，全局坐标系下
    const auto& qface_points = fe_elem_face->get_xyz();

    // 单元的residual和jacobian
    DenseMatrix<Number> Ke;
    DenseVector<Number> Fe;

    // 左变量对左单元residual的jacobian
    DenseMatrix<Number> Kee;
    // 右变量对左单元residual的jacobian    
    DenseMatrix<Number> Kne;
    // 左变量对右单元residual的jacobian    
    DenseMatrix<Number> Ken;
    // 右变量对右单元residual的jacobian    
    DenseMatrix<Number> Knn;

	DenseVector<Number> Fee;
	DenseVector<Number> Fen;

    // 当前单元的自由度索引
    std::vector<dof_id_type> dof_indices;

    // 罚值
    // TODO: 把lamda提出到配置文件中
    Real lamda = 1E+0;

    // 单元平均值
    // TODO: 引入真实的单元平均值
    // Real average = 1;

	NumericVector<Real>& solution  = es.parameters.get<FEProblemBase*>("_fe_problem_base")->getNonlinearSystemBase().solution();
	mooseInfo(solution.local_size());
	for (size_t i = 0; i < solution.size(); i++)
	{
		//std::cout << solution(i) << std::endl;
	}

	std::vector<DenseVector<Number>> phi_average;
	for (const auto& elem : mesh.active_local_element_ptr_range())
	{
		auto id = solution(elem->id());
		// 更新elem的索引
		dof_map.dof_indices(elem, dof_indices);
		auto n_dofs = dof_indices.size();

		// 更新elem的fe
		fe->reinit(elem);

		DenseVector<Number> loc_phi_average;
		loc_phi_average.resize(n_dofs);
		Real elem_volume = 0;

		for (unsigned int _qp = 0; _qp < qrule.n_points(); _qp++)
		{
			elem_volume += JxW[_qp];
			for (unsigned int _i = 0; _i < n_dofs; _i++)  // 基函数循环
			{
				loc_phi_average(_i) += lamda * JxW[_qp] * phi[_i][_qp];
			}
		}

		for (unsigned int _i = 0; _i < n_dofs; _i++)  // 基函数循环
		{
			loc_phi_average(_i) /= elem_volume;
			std::cout << loc_phi_average(_i) << std::endl;
		}
		phi_average.push_back(loc_phi_average);
	}

    // 单元循环
    for(const auto& elem : mesh.active_local_element_ptr_range())
    {
		auto id = elem->id();
        auto average = solution(elem->id());
        // 更新elem的索引
        dof_map.dof_indices(elem, dof_indices);
        auto n_dofs = dof_indices.size();

        // 更新elem的fe
        fe->reinit(elem);

        // 开辟residual和jacobian空间
        Ke.resize (n_dofs, n_dofs);
        Fe.resize (n_dofs);

        // 面循环
        for (auto side : elem->side_index_range())
         {
            // elem不存在neighbor，为边界面
            // TODO: VFV重构暂不引入边界条件，如何定义边界条件？
            if (elem->neighbor_ptr(side) == nullptr)
             {
                // 更新Face上的fe
                fe_elem_face->reinit(elem, side);

                // 创建cell的face
                std::unique_ptr<const Elem> elem_side (elem->build_side_ptr(side));

                for (unsigned int qp=0; qp<qface.n_points(); qp++)
                  {
                    for (unsigned int i=0; i<n_dofs; i++)
                     {
                        // 边界的resiudal
                        //Fe(i) += 0;
                        for (unsigned int j=0; j<n_dofs; j++)
                         {
                            // 边界的jacobian
                           // Ke(i,j) += 0;
                         }
                     }
                 }
             }

            //内部面
            else
              {
                // 邻单元
                const Elem * neighbor = elem->neighbor_ptr(side);

                // 获取单元编号
                const dof_id_type elem_id = elem->id(), neighbor_id = neighbor->id();

                // 判断当前面是否参与计算
                // elem_id < neighbor_id 为了避免面被计算两次
                if ((neighbor->active() && (neighbor->level() == elem->level()) && (elem_id < neighbor_id)) ||
                    (neighbor->level() < elem->level()))
                  {
                    // 创建当前Face
                    std::unique_ptr<const Elem> elem_side (elem->build_side_ptr(side));

                    // 更新单元在当前面上的fe，如基函数
                    fe_elem_face->reinit(elem, side);

                    // 获取面的积分点，全局坐标系下的
                    std::vector<Point > qface_point = fe_elem_face->get_xyz();

                    // 定义邻单元在当前面积分点坐标，局部坐标系下的
                    std::vector<Point> qface_neighbor_point;

                    // 计算积分点坐标，全局坐标系-->局部坐标系下的反变换
                    FEInterface::inverse_map (elem->dim(),
                                                fe->get_fe_type(),
                                                neighbor,
                                                qface_point,
                                                qface_neighbor_point);

                    // 更新邻单元在当前面`qface_neighbor_point`上的fe，如基函数
                    fe_neighbor_face->reinit(neighbor, &qface_neighbor_point);

                    // 获取邻单元的自由度索引
                    std::vector<dof_id_type> neighbor_dof_indices;
                    dof_map.dof_indices (neighbor, neighbor_dof_indices);
                    const unsigned int n_neighbor_dofs = cast_int<unsigned int>(neighbor_dof_indices.size());

                    // 开辟Face的jacobian
                    Kne.resize (n_neighbor_dofs, n_dofs);
                    Ken.resize (n_dofs, n_neighbor_dofs);
                    Kee.resize (n_dofs, n_dofs);
                    Knn.resize (n_neighbor_dofs, n_neighbor_dofs);

					Fee.resize(n_dofs);
					Fen.resize(n_neighbor_dofs);
                    // face上的Gauss积分循环
                    for (unsigned int qp=0; qp<qface.n_points(); qp++)
                    {
						for (unsigned int i = 0; i < n_dofs; i++)   
						{
							Fee(i) += 0; //lamda * JxW_face[qp]*(solution(elem_id) - solution(neighbor_id))*(phi_face[i][qp]-phi_average[elem_id](i));
						}

						for (unsigned int i = 0; i < n_neighbor_dofs; i++)
						{
							Fen(i) += 0; // lamda* JxW_face[qp] * (solution(elem_id) - solution(neighbor_id))* (phi_neighbor_face[i][qp] - phi_average[neighbor_id](i));
						}

                        for (unsigned int i=0; i<n_dofs; i++)  // 基函数循环
                        {
							for (unsigned int j=0; j<n_dofs; j++)  // 测试函数循环
                            {
								Kee(i,j) += JxW_face[qp]*(phi_face[i][qp]-phi_average[elem_id](i))*(phi_face[j][qp]);  //phi(+) * test(+)
                            }
                        }

                        for (unsigned int i=0; i<n_neighbor_dofs; i++)
                        {
							for (unsigned int j=0; j<n_neighbor_dofs; j++)
                            {
								Knn(i,j) += JxW_face[qp]*(phi_neighbor_face[i][qp]- phi_average[neighbor_id](i))*(phi_neighbor_face[j][qp]);  // -phi(-) * test(-)
                            }
                        }

                        for (unsigned int i=0; i< n_dofs; i++)
                          {
                            for (unsigned int j=0; j< n_neighbor_dofs; j++)
                              {
                                Ken(i,j) -= JxW_face[qp]*(phi_face[i][qp]- phi_average[elem_id](i))*(phi_neighbor_face[j][qp]);  // -phi(+) * phi(-)
                              }
                          }

                        for (unsigned int i=0; i<n_dofs; i++)
                         {
                            for (unsigned int j=0; j<n_neighbor_dofs; j++)
                             {
                                Kne(i,j) -= JxW_face[qp]* (phi_neighbor_face[i][qp] - phi_average[neighbor_id](i))*(phi_face[j][qp]); // phi(+) * phi(-)
                             }
                         }
                     }

                    // 将Kne，Ken，Kee，Knn加入到LineaImplictSystem的稀疏矩阵中
                    _sys.matrix->add_matrix(Kne, neighbor_dof_indices, dof_indices);
                    _sys.matrix->add_matrix(Ken, dof_indices, neighbor_dof_indices);
                    _sys.matrix->add_matrix(Kee, dof_indices);
                    _sys.matrix->add_matrix(Knn, neighbor_dof_indices);
					_sys.rhs->add_vector(Fee, dof_indices);
					_sys.rhs->add_vector(Fen, neighbor_dof_indices);

					//std::cout << Kne << std::endl;
					//std::cout << Ken << std::endl;
					//std::cout << Kee << std::endl;
					//std::cout << Knn << std::endl;

					//mooseError("Here");
                 }
             }
         }
    }     
    auto end = steady_clock::now();
    duration<double> time_span = end-start;
    Moose::out << "组装消耗的时间：" << time_span.count() << "\n"; 
}

// ReconstructionSystem ////////

ReconstructionSystem::ReconstructionSystem(FEProblemBase & subproblem, const std::string & name) :
	SystemBase(subproblem, name, Moose::VAR_AUXILIARY),
	PerfGraphInterface(subproblem.getMooseApp().perfGraph(), "ReconstructionSystem"),
    _fe_problem_base(subproblem),
	//_vfv_problem(subproblem),
    _mesh(subproblem.mesh()),
    _eq(subproblem.es()),
    _sys(subproblem.es().add_system<LinearImplicitSystem>(name)),
	_assembly(_fe_problem_base.assembly(0)),
	_current_solution(nullptr),
	_serialized_solution(*NumericVector<Number>::build(_fe_problem_base.comm()).release()),
	_solution_previous_nl(nullptr),
	_u_dot(nullptr),
	_u_dotdot(nullptr),
	_u_dot_old(nullptr),
	_u_dotdot_old(nullptr),
	_need_serialized_solution(false)
{
	_elem_vars.resize(libMesh::n_threads());

}

ReconstructionSystem::~ReconstructionSystem() { delete& _serialized_solution; }

void ReconstructionSystem::init()
{
	// 调用System::update(), 同步processors的数据
	this->update();

  // 什么都不做，占位
  //_sys.add_variable("u_reconstruction", _fe_type);
  //_sys.get_linear_solver()->set_solver_type(_solver_type);
  //_sys.get_linear_solver()->set_preconditioner_type(_preconditioner_type);
}

void ReconstructionSystem::initialSetup()
{
	mooseInfo("ReconstructionSystem::initialSetup");
	this->solve();
}

void ReconstructionSystem::timestepSetup()
{
}

void ReconstructionSystem::subdomainSetup()
{
}

void ReconstructionSystem::residualSetup()
{
}

void ReconstructionSystem::jacobianSetup()
{
}

void ReconstructionSystem::updateActive(THREAD_ID tid)
{
}

void ReconstructionSystem::reinitElem(const Elem* elem, THREAD_ID tid)
{
	mooseInfo("ReconstructionSystem::reinitElem");
	const std::vector<MooseVariableFEBase*>& vars = _vars[tid].fieldVariables();
	MooseVariableFEBase* var = vars[0];
	var->prepare();
	var->computeElemValues();

}

void ReconstructionSystem::reinitElemFace(const Elem* elem, unsigned int side, BoundaryID bnd_id, THREAD_ID tid)
{
}

NumericVector<Number>&
ReconstructionSystem::serializedSolution()
{
	if (!_serialized_solution.initialized())
		_serialized_solution.init(_sys.n_dofs(), false, SERIAL);

	_need_serialized_solution = true;
	return _serialized_solution;
}

void ReconstructionSystem::augmentSparsity(SparsityPattern::Graph&, std::vector<dof_id_type>&, std::vector<dof_id_type>&)
{
}

void
ReconstructionSystem::serializeSolution()
{
	if (_need_serialized_solution &&
		_sys.n_dofs() > 0) // libMesh does not like serializing of empty vectors
	{
		if (!_serialized_solution.initialized() || _serialized_solution.size() != _sys.n_dofs())
		{
			_serialized_solution.clear();
			_serialized_solution.init(_sys.n_dofs(), false, SERIAL);
		}

		solution().localize(_serialized_solution);
	}
}

void ReconstructionSystem::addVariable(const std::string & var_type, const std::string & name, InputParameters & parameters)
{
	parameters.set<MooseEnum>("order") = "FIRST";
	SystemBase::addVariable(var_type, name+"_rs", parameters);
}

//void ReconstructionSystem::addVariable(const std::string& name, FEType type)
//{
//	setFEType(type);
//	_sys.add_variable(name, type);
//}

void ReconstructionSystem::solve()
{

	Moose::out << "init \n" ;
	_sys.attach_assemble_function(assemble_reconstrunction);
    _fe_problem_base.es().reinit();

    Moose::out << _sys.get_linear_solver()->solver_type() << '\n';
    auto start = steady_clock::now();
    _sys.solve();
    
    auto end = steady_clock::now();
    duration<double> time_span = end-start;
    Moose::out << "求解消耗时间：" << time_span.count() << "\n";
    
    Moose::out  << _eq.parameters.get<Real>("linear solver tolerance") << '\n';

    const unsigned int n_linear_iterations = _sys.n_linear_iterations();
    const Real final_linear_residual = _sys.final_linear_residual();
    
    Moose::out << "线性迭代步数: " <<  n_linear_iterations << '\n';
    Moose::out << "残值: " <<  final_linear_residual << '\n';

    ExodusII_IO (_eq.get_mesh()).write_equation_systems ("out.e", _eq);
}